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£\j ■ We consider analytical and numerical solution of NMR relaxation under the condition of surface 

relaxation in an equilateral triangular geometry. We present an analytical expression for the Green's 
function in this geometry. We calculate the transverse magnetic relaxation without magnetic gradi- 
ents present, single-phase, both analytically and numerically. There is a very good match between 
the analytical and numerical results. We also show that the magnetic signal from an equilateral 
triangular geometry is qualitatively different from the known solution: plate, cylinder and sphere, 
in the case of a nonuniform initial magnetization. Non uniform magnetization close to the sharp 
corners makes the magnetic signal very fast multi exponential. This type of initial configuration 
fits qualitatively with the experimental results by Song et al. [Jj. It should also be noted that the 
solution presented here can be used to describe absorption of a chemical substance in an equilateral 
^ ■ triangular geometry (for a stationary fluid). 
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I. INTRODUCTION 
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£\\ • In the oil industry there has been a great interest in using NMR as a tool for improved reservoir characterizing. 
t-H ' NMR can be used as an in situ tool for measuring oil and brine content in saturated porous rocksjlllllll. The 
distribution of oil and water in the porous rock is also very important. If a phase is in contact with pore wall, it 
will experience enhanced relaxation]^ Q- One way of observing this is by studying how the oil signal changes when 
the core is aged, see H, |t| and the references therein. The distribution of oil and water in a porous rock is mainly 
determined by the geometry of the pore space and the chemical composition of oil, brine and rock. The NMR signal 
from a porous rock will then be a strong function of the specific pore geometry and the wetting status of the rock. 
Experiments on porous rocks are well suited for showing an effect of enhanced surface relaxation due to wettability 
change. However, by study simpler systems, more information can be gained in order to understand how chemical 
I • properties will affect the magnetic signal. 

By assuming a specific geometry for the porespace more information from the NMR signal can be obtained. Cylin- 
drical and spherical geometries can be used when wettability is not a topic. However most reservoirs are mixed wet 
and one need geometries containing sharp corners. This realization is the starting point for this work. One of the 
• • . simplest geometries which allows for more than one phase to form a stable configuration is an equilateral triangular 
geometry. A great deal of attention have been devoted to triangles in order to understand the multi phase behavior 
, of porous rockspS [H Il2| . Triangles are the basic building blocks in pore network models[l3l| and also bundle of 
5_i 1 equilateral triangular tubes models have been studied extensively [fH llq . Attempt of using a bundle of equilateral 
triangular tubes in interpreting NMR signal from porous rocks have also been done [HI- 

As always with NMR experiments the key challenge is to get a proper interpretation of the magnetic signal. The 
crucial step in order to gain precise knowledge from NMR signal is to have a good theoretical description of the 
magnetic signal. This work will only be concerned with magnetic signal from single phases, but it forms a necessary 
basis in developing numerical code for magnetic signal from triangles containing multi phases. The numerical solution 
presented in this work can be extended to account for more than one phase inside the rock and analytical single 
phase results is necessary in order to calibrate the numerical algorithm. The multi phase numerical code can then 
be used to study how wettability will influence the magnetic signal in a well defined geometry. Experiments on glass 
triangular tubes are in progress at our group and the analytical solution presented here will be used to interpret the 
experimental result, this will be considered in a different work. 

The higher modes presented here are usually neglected, but have been shown to give valuable information of pore 
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sizes [H Il7|. In order to fully characterize a geometry one need information from higher modes. The lowest mode is 
dependent on the surface relaxivity and proportional to the internal surface area divided by the total volume of the 
pore (S/V), the higher modes are independent of the surface relaxivity and proportional to (S/V) 2 . A clever way of 
enhancing the weight of the higher modes is by creating a non uniform initial magnetization, as described by p|. 

The outline of the papers is as follows: in section [H] we present the theory for magnetic relaxation in a confined 
space. In section Hill we calculate the Green's function for an equilateral triangular geometry and in section ITvl we 
calculate the magnetization for different initial configurations. Then we present details of the random walk solution 
in section Ivl and finally conclusion and discussion. 

II. THEORY 

The equilateral triangular geometry is a true two dimensional system, contrary to plate, cylinder and sphere 
geometry, first considered by Brownstein and Tarr0. This problem is much harder to solve, the technical reason for 
this is that it is not possible to choose a coordinate system where the axes are parallel to the sides of the triangle. 
Fortunately the calculation greatly simplifies because of a series of recent papers by Mccartin, which solves the diffusion 
equation in an equilateral triangular geometry [IH ITii l2(il]. The magnetization as a function of time is given by the 
following equations @, 0: 

2MM = i ,v^(M)-^. a) 

T 2f) is the bulk relaxation and D the diffusion constant. At the surface, S, we have the following boundary condition: 

Dn-VM(r,t)+pM(r,t)| r=E = 0, (2) 

where n is the outward normal and p the surface relaxivity. Equation (JJJ (with I/T2& = 0) and J5J) are identical 
to the equations governing absorption of a chemical substance on a surface, in the case of a stationary fluid. The 
magnetization of the sample is : 

M(t) = [ drM(r,t) , (3) 

J A 

where the integral is taken over the triangular domain. By introducing the Green's function we can write the 
magnetization as: 

M(r, t) = e- t/T2b [ dr' p(r')G(r'\r;t) , (4) 

J A 

where p(r') is the initial spin density. G(r'|r;£) is the Green's function or the propagator. It is defined as the 
probability for a particle at position r' at time to diffuse to point r during a time t. The propagator satisfies the 
diffusion equation at the interior of the pore space: 

dG{ f t ' ]t) -DV 2 G(v\r>;t) = 0; and 

G(r|rV)|t=o=5(r-r'), (5) 
where D is the diffusion constant. The boundary condition at the surface S : 

Dh- VG(r|r';i) + pG(r|r';i)|r=s =0. (6) 



III. THE GREEN'S FUNCTION IN AN EQUILATERAL TRIANGULAR GEOMETRY 

In this section we will discuss relevant eigenf unctions for situations where there is no gradients present, for the 
general case see Appendix^ By using the standard eigenfunction expansion of the propagator: 

00 

G(r\r';t)=J2Mr)Mr')e- t/T ', (7) 

4=0 



FIG. 1: Equilateral triangle of side length a and inscribed radius R 



where {4>i\ are an orthonormal set of eigenfunctions with corresponding eigenvalues T^. From equation J^} and © it 
then follows: 



L»V 2 ^(r) = -— <Pi(r) and 
^ i 

Dh- V0 i (r)+p^.(r)| r=s = O. 



(8) 



This equation needs to be solved in an equilateral triangle, shown in Figure ^ The orthonormal set of eigen- 
functions can be divided in a set of eigenfunctions symmetric and antisymmetric about the line x = a/2, i.e. 
{(/,} = cf) a } = {T s /N s , T a /N a }. {T s ,T a } is the set of orthogonal eigenfunctions and N s > a the normalization. 
When the magnetization is given by equation the only relevant modes are the symmetric diagonal modes (for the 
full solution see Appendix 0) : 
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and the eigenvalues 
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\x and 62 are determined from the boundary condition JSJ and it give rise to the following transcendental equation: 
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tan7r/i = , jj, e< 1, 1 + *l 

27T/i 
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7 : 



pi? 



2tt/x ' D 

A reduced from of the propagator, when the magnetic signal is given by equation Q is then: 



G(x,y\x',y';t) = ^ 



Tl{x,y)Tfi{x'> y ') „_ m 



i=0 



(11) 



(12) 



This equation can now be used to calculate the magnetic signal from an equilateral triangular pore given an initial 
spin density. 
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TABLE I: Different initial spin densities 



Type 




uniform 


1/(3a/3#0 


center 


S(x- V3R) 5(y-R) 


side 


5(x- V3R) 5{y- R/32) 


corner 


8(x- V3R) 6{y- 7 R/32) 



Triangle - uniform initial spin density 



Triangle - center initial spin density 



M n.f 
1 n.r 











1 1 1 
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FIG. 2: Upper-Left: the first four intensities for uniform initial spin density (Jq — > 0.61, Ji — > 0.15, J2 — ► 0.07, J3 — ► 0.04). 
Upper-right: the first four modes for all spins initially at the center of the triangle (Jo — > 1.65, I\ — * —0.83, I2 — > 0, J3 — > 0.41). 
Lower-left: the first four modes for all spins initially at the side of the triangle (Jo — > 0.94, Ji — » —0.13, I2 — > 0.54, J3 — > —0.32). 
Lower-right: the first four modes for all spins initially at the corner of the triangle (Jo — ► 0.04, Ji — > 0.13, I2 — > 0.24, J3 — > 0.32). 



IV. MAGNETIC SIGNAL FROM AN EQUILATERAL TRIANGLE 



From equation ©, J3J and lfP2)l we find : 



M(t) =e- t/T2b J2 / drdrV(rO0■(r)^(r>- ^ / T ' 
e -*/T„ 



2 ^iVl(^) 2 



dr'p(r')i;i(r') (cos[c5 2 ] - cos[5 2 - 2//tt] + 2^sin[<S 2 ]) 



(13) 



We need to know the initial excited spin density. Usually it is assumed that the initial spin density is uniform as 
in the work by Brownstein and TarrQ. However, recently Song et al.Q] have shown during a series of papers that 
it is possible to create a nonuniform initial magnetization. The physical explanation for this is due to susceptibility 
differences between water and rockQ. There are then four interesting cases, shown in tabled For the side and corner 
densities we have chosen points close to the surface, which lies at one of the symmetry line of the triangle. Then we 
only need to consider one corner (side) and the result will be equivalent for the other corners (sides). 
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(14) 
(15) 

(16) 
(17) 



where is given in equation 0. 

We have compared the analytical solution for the equilateral triangle with the solution for plate, cylinder and 
sphere, with different initial spin densities. Initial spin density close to a corner clearly does not exists for those 
geometries. We do not show the result for plate, cylinder and sphere here, but they behave qualitatively in the 
same manner. When approaching the slow diffusion limit (SDL, see next subsection) the weight of the higher modes 
increases compared to the uniform case, but the lowest mode always gives the highest contribution to the magnetic 
signal, see Figure [2] From Figure [3 it is evident that initial spin density close to a corner is qualitatively different 
from uniform, center and side. Sharp corners and initial spin density close to the corners make the magnetic signal 
multi exponential very fast, even for 7 = 1. For 7 = 4 the lowest mode is the least dominant one. 

In order for the lowest mode to be dominant the spin density for all times must be more or less uniform. In the 
non uniform cases, the spins lose coherence due to surface effects before the magnetization is uniform. This has the 
effect that non uniform initial spin density in general makes the higher modes more dominant. In order to suppress 
the lowest mode, the surface effects must be felt before the spins manage to diffuse over the triangular domain to 
give a uniform magnetization. Clearly, corners will be very efficient of suppressing the lowest mode. This also give a 
natural explanation as to why this effect is not so prominent for plate, cylinder and sphere. 



A. Fast Diffusion and Slow Diffusion Limit 



There are two different time scales of interest. The relaxation time which is dependent on the surface relaxivity, 
r p ~ R/p and the diffusion time which is dependent on the diffusion constant td ~ R 2 / D. If m « t p we are in the 
fast diffusion limit (FDL) an if f r p << td then we are in the slow diffusion limit (SDL). In these two limiting cases 
the eigenvalues and the analytical expressions for the magnetic signal simplify considerably. In the FDL the spins 
traverse the triangular domain many times before they relax and the magnetic decay is dominated by one mode. By 
simply replacing tan with its argument in equation (|1 1|) . we regain the famous result by Brownstein and Tarr0, by 
direct calculation: 

12 5 

5E-"s->7- (18) 

We write the magnetic signal from the equilateral triangle as M(t) = In exp(—t/Tu), the coefficients and eigen- 
values are summarized in table |n] and IIIII As seen from table IIIII the FDL is dominated by one mode and hence a 
single decay time. In Figure [3] we have plotted the eigenvalues as a function of 7. 

Figure and |21 show that the eigenvalues deviate faster from the FDL compared to the intensities created with an 
uniform initial spin density. This means that the magnetic decay is still mono exponential, but the simple relation 
T^ 1 = pS/V is violated. 



V. NUMERICAL ALGORITHM - RANDOM WALK 



In this section we will show explicitly that a random walk algorithm]^ gives rise to the diffusion equation QJ, 
with the appropriate boundary condition given in equation |2"|l. when the number of grid points approaches infinity. 



6 



TABLE II: Characteristic decay times as a function of surface relaxtivity , pore radius R and diffusion coefficient D in an 
equilateral triangle. Note that the lowest mode is independent of the diffusion constant in FDL, while the higher modes are 
independent of the surface relaxtivity 





Fast Diffusion Limit 
TD » T p (pR/D « 1) 


Slow Diffusion Limit 

T D « T p {pR/D » 1) 


(22)00 


R 

^ 


9R 2 
4Dir' 2 


(T 2 )u 


9R 2 
iDTr'-'i 2 


9R A 



TABLE III: The intensities in the fast and slow diffusion limit. Note that the sum of the intensities add up to one as they 
should, e.g. J2Zo + if = ^ '/6 and J2Zo sin [2(1 + »)t/3]/(1 +i) = tt/6. In FDL (r p >> td) hi = 5 ifi for all initial spin 
densities. 





Slow Diffusion Limit 

TD » T p (pR/D » I) 


In uniform 


6 

(l + i) 2 n 2 


In center 


6sinl2(l+i)7r/3J 
(l+»)ir 


la side 


4siii|(i + i)7r/8|+2sin|V(l + i)7r/4| 


(1+Oir 


la corner 


4sin|7(l-))7T/S +2 sin[(l + i)7r/4J 


(l + i)7T 



The random walk algorithm presented here will serve two purposes: (1) it will give a confirmation on the analytical 
results and (2) it can be extended to account for any number of phases inside a triangle. More details on scaling and 
convergence in random walk simulations are given in appendix FBI 

For random walk simulation in an equilateral triangle we choose a hexagonal grid, see Figure 0] Random walkers 
are placed at random at a lattice point in the equilateral triangle. The number of random walkers at an interior point 
r ( see Figure when the clock advance one step r is then given by : 

M(r,t + r) - M(r,t) = 

- \M(r + ei, t) - M(r, t) + M (r - el, t) - M (r, t) 
6 L 

+ M(r + ej, t) - M (r, t) + M(r - ej, t) - M (r, t) 
+ M(r + efc, t) - M(r, t) + M (r - efe, t) - M (r, t) 

-nM(r,t) (19) 



The probability for a random walker to take a step in each direction is 1/6, the probability for a random walker to 
die during the time step r is k. Dividing by r, we find: 

M(r,t + r) - M(r,t) _ 



e 

g7 



M (r + ei, t) - 2M (r, t) + M(r - ei, t) 



M (r + ej.t) - 2M (r, t) +M(r- ej, t) 



M (r + ek, t) - 2M (r, t) +M(r- ek, t) 



-M(r,t) 



(20) 



The lattice spacing is given by e. Taking the limit r 



dM 
~dt 



e 

6r 



d 2 M 
di 2 



+ 0, e- 



0, K - 

d 2 M 

dk 2 



- -M(r,t). 

T 



(21) 



7 




~~ i — i — i — I — r 







n 10 20 30 40 50 60 70 80 90 100 
7 = pR/D 



To/Ifdl - 
r 2 /T FDl - 










/ / k 








* ' "1 


t 

/ 



/ a 
,' .<*' 

' 









0.01 0.1 1 10 100 

7 = pB/D 



r 2 /r FDL 



FIG. 3: Left: The eigenvalues as a function of 7, Tfdl = r/(2p). Right: same as left but log- log plot. For 7 = 0.1 there is a 
10% deviation and at 7 = 0.4 there is a 20% deviation from the FDL. 




FIG. 4: Grid for random walk simulation, lattice spacing e and side length a, unit outward normal vectors 61,2,3, 



Changing from lattice coordinates to Cartesian coordinates: 



di 2 + df 



dk 2 



3 (&_ 

2 \dx 2 



dy 2 



(22) 



we finally arrive at 



— = D\7 2 M - — , where 
at T 2b 



° = Tr 



and 



1 



(23) 
(24) 



This is the same equation as Q. At a lattice point at a boundary surface normal to the n3-vector (in Figure 0J), the 
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equation becomes: 

M(r,t + r) - M(r,t) 



c- 



2 f_ 
3 e 4r 



M (r + ek, t) - 2M(r, t) + M(r - ek, t) 



e 2 



2 e 2 M(r - ei,t) - M(r,t) 
+ 3 47 e 

2 e 2 M(r - ej,t) - M(r,t) 
+ 3~ir e 

-~M(r,t) 

3 T 

-e-M(r,t). (25) 
r 

We have assumed that the walkers have a probability ( of being killed when stepping from an interior point to the 
wall. If a walker is not killed it is assumed to return to the interior point in the same time step. In the limit r — ► 0, 
e — > 0, the fractions involving M are recognized to approach dM/dt, d 2 M/de 2 , —dM/di and —dM/dj, respectively. 
Using the following relation: 

| + | = v^n a .V (26) 

and as the LHS and the first term on the RHS in equation (|25H arc of higher order in e, they can be neglected compared 
to the others, hence : 

Dn 3 • VM + pM = and p = — ^ . (27) 

2\/3 T 

By symmetry the other boundary sides give the same answer and we then regain equation J3J). 
Using the fact that number of lattice points from corner to corner along an edge is 

N = l + a/e, (28) 

from equation l|24|l and (|27|l we find the following important relation: 

l= E § = l(N-l)C (29) 

For a given 7 the probability for a random walker to die when it hits the wall can be calculated from eauation (|29|l . 
Consider the ratio of bulk lattice points to surface lattice points we can define a new parameter 6, which holds the 
information of bulk relaxation T 2 6 : 

N 1 + 1/N k_ 1 a N 2 (N + 1) 
T(l-l/iV) 2 C ~ 4^ P T 2b (N — l) 3 

(30) 



1 

1 4D 



The magnetic signal can then be calculated numerically by placing a number of random walkers inside the equilateral 
triangular domain. For uniform initially spin density the walkers are placed at random and for the delta densities all 
the walkers start in a point of the triangle. At each time step the walkers take one step of length e and dies with a 
probability n. If the walkers hit the wall they die with probability £. The magnetic signal will then be proportional 
to the number of walkers alive at each time step r. The lattice spacing and time is related to the physical length and 
time by using equations H24|l. 

In Figure and we have compared numerical and analytical results. There is clearly a very good match. The 
largest deviation between numerical and analytical results are for long times. This is natural as the number of walkers 
is low and the statistic is poorer. Note that if one extrapolate the straight line for high 7 values in Figure [5] (left) and 
El (left), it crosses the y-axis at 6/ir 2 and 6/77, which is expected from Table ITTT1 We have also made some comparison 
with 9 ^ 0, i.e. with bulk relaxation (see equation |j3U|)). In Figure [7| we have plotted results for 6 = 0.1 and 1. For 
a = 100p.ni, D — 2500/xm 2 /s and 7 = 1, 6 = 0.1, 1 corresponds to a bulk relaxation of 1.7 and 0.17 s respectively. 
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FIG. 5: Left: Comparison of random walk simulation (points) and analytical solution (line) for 7 = 4 and 100. TV = 501 and 
10 6 random walkers. The initial spin density is uniform. Right: Same as left but log scale on the x-axis instead of the y-axis. 




FIG. 6: Left: Comparison of random walk simulation (points) and analytical solution (line) for 7 = 0.1, 1 and 10. The initial 
spin density is the center density. TV = 501 and 10 6 random walkers. Right: Same as left but log scale on the x-axis instead of 
the y-axis. 



A. Scaling in random walk simulations 

Eliminating e between equation 1|24[) and equation l|28|) gives a relation between the intrinsic time scale for diffusion 
and the time step length in a simulation with a given TV value: 

r D ^r(N-lf = ^ (31) 

If, for given 7 and 9 values, the average magnetization is plotted not against the actual time t or the integer time step 
counter t/r but in units of the intrinsic time scale, viz., a plot against 

t _ t 
Td ~ t(TV- l) 2 

then the results for various {TV, £} consistent with the given 7 value should be expected to plot on top of each other, 
provided inaccuracies due to finite TV values are unimportant. This follows since t/r oc (TV — l) 2 , where TV is not a 
physical parameter but an artificial one determined by the simulation conditions. Technically, there would thus be a 
scaling property in the limit TV — > 00, in that the magnetization should depend only on a certain algebraic combination 
of simulation parameters. 

The predicted scaling has been checked numerically for several 7 values, using 9 — (no bulk relaxation), and has 
indeed been found to hold. (The origins and magnitude of deviations from exact scaling is discussed in more detail 
in Appendix[51) Its practical value is that for a given simulation, one may choose that set of values for {TV, £} which 
gives an optimally acceptable combination of accuracy and computing time requirements. This scaling is also very 
useful when comparing numerical and analytical results, once the results have been normalized by t — > t/rn only one 
parameter, 7, has to be given. Pore size a, diffusion constant D and surface relaxtivity p are unimportant as long as 
the combination 7 = pR/D — pa/ stay unchanged. 

This scaling is also consistent with the one found by Valfouskaya et al. |24ll25l.l2i| . They show that for zero surface 
relaxtivity an universal curve exists for a porous media. It is universal in the sense that different saturations give 
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FIG. 7: Left: Comparison of random walk simulation (points) and analytical solution (line) for 9 — 0.1, 1, 7 = 1. Uniform 
initially spin density, N = 201 and 10 6 random walkers. Right: Same as left but log scale on the x-axis instead of the y-axis. 
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FIG. 8: Magnetic signal for uniform and corner initial spin density (7 = 1). The points are random walk calculations and lines 
analytical. 



the same curve. Considering magnetic signal from a single pore, the curve by Valfouskaya et al. reduce to the 
effective diffusion constant normalized by the molecular diffusion constant (D), D e g(Ti) / D = (r(T]) 2 ) /(4iD), where 
Ti = 4/(9^/W)\nji/l s . l s is a typical length scale of the porous media, in our case it is ^ s = V/S = a/(4"\/3) and 
T[ = 2/(9y / 7r)y / t/rD, where td is defined in equation lpTT|) . We have showed, in the beginning of this subsection, that 
once 7 is fixed (and zero bulk relaxtivity) all magnetization curves, independent of the initial spin density, will plot 
on top of each other when the physical time has been rescaled according to t — > t/T]j. This is exactly the type of 
scaling used bv plEll^. 



VI. DISCUSSION 



In pore network modeling literature, the introduction of sharp corners has shown to be very fruitful in order to 
understand multi phase behavior. In this paper we have considered magnetic signal from one of the simplest geometries 
containing sharp corners, the equilateral triangle. This solution is similar to the one dimensional geometries: plate, 
cylinder and sphere in the case of uniform initial magnetization. When the initial magnetization is non uniform 
the triangle behaves qualitatively different. In Q it was shown experimentally how to generate a non uniform initial 
magnetization, due to diffusion in internal field (DDIF). A reference signal with uniform initial magnetization was used 
and a signal due to a non uniform initial magnetization was created. It was seen that the non uniform magnetization 
was highly multi exponential. With the triangular geometry this result can be explained by calculating the magnetic 
signal with uniform initial magnetization and magnetization concentrated close to the corners, shown in Figure |S] 
Compared with plate, cylinder and sphere only the triangle with non uniform initial magnetization close to the corners 
had a magnetic signal which behaved multi exponential for 7 = 1. 

In order to fully explain the experimental results generated by the DDIF technique, effect of pore connectivity |27| 
and internal gradients must also be considered. However, the results shown in Figure [5] are very striking and they 
seem to capture some of the essence in the experimental results^]. 

The numerical results presented here can be extended to account for more than one phase inside the triangle. 
These results will then be used to study the influence of wettability on the magnetic signal and will form a basis when 
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constructing inversion integrals used to interpret NMR-logs. 



VII. CONCLUSION 

We have presented an analytical solution for the magnetic signal from an equilateral triangular pore, with surface 
relaxation. To our knowledge this solution has not been presented before. This solution will be used in theoretical 
studies of single- and two-phase NMR response from equilateral triangles, which can be used as basic building blocks 
for pore scale models. This solution is also very important when interpreting experiments done at equilateral triangular 
tubes, which is in progress at our group. 

We have studied how the magnetic signal changes for different initial conditions. For non uniform initial magneti- 
zation the equilateral triangular geometry behaves qualitatively different than plate, cylinder and sphere due to the 
sharp corners. 

In a forthcoming paper we will also present results for the pulsed field gradient spin echo (PFGSE) sequence psl I29I 
IscLlSlj . the single-phase result can be calculated from the Green's function presented in this paper. Two-phase results 
for any value of 7 must in general be solved by a random walk algorithm, except for the limit 7 — > where analytical 
results can be found. The analytical result will then serve to calibrate the numerical random walk algorithms for 
two-phase. 
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APPENDIX A: THE FULL GREEN'S FUNCTION FOR AN EQUILATERAL TRIANGLE 



The Green's function for the equilateral triangle is given by the following equation. 

00 

G(x, y\x',y'; t) = £ <^(x, yW ht (x' , y')^'^ 



OO OO 



■£ E [C■(^ J y)< J ■(^^yO + C•^^)< J '( a: ^y')]e- ^/T - 



i=0 i=i+l 

rpS,a 

k s,a _ i,j 



C?° = ttf^ and N *<? = / dxdy T *f( x > v ^ T if( x > y) ■ 

i,j J A 



(Al) 
(A2) 



The complete set of orthogonal eigenfunctions for the equilateral triangular triangle are|20j: 



T ii(x,y) = cos 



ttA 
3R 



(3R-y)-S 1 



Vo7r(/x — v) 



9R 



(x - V3R) 



+ cos 



cos 



VStt(v - A) 



T?Ax,y) = cos 



ttA 
3^ 



(3R - y) - <5i 



9R 



sm 



[x - VZR) 



cos 



^(3R-y)-S 3 



cos 



v37r(A — /i) 



(A3) 



9R 



(x - V3R) 



V37r(/i — v) 



9R 



(x - V3R) 



(A4) 



-I- cos 



^(3i?-y)-<5 2 



sm 



V3tt(i< - A) 



9R 



(x - V3i?) 



-I- cos 



^(3R-y)-6 3 



sm 



y/3n(\ — j-i) 



9R 



(x - V3R) 



R = a/(2\/3) is the radius of the inscribed circle in an equilateral triangle of side a. The index s, a of the eigenfunctions 
refers to symmetric, antisymmetric about the line x = a/2, respectively. The complete set of orthogonal eigenfunctions 
is {T?j(i > j); T^M > j)}, for further details on this point seej^- The constants fi, v, A, 5\, 62, S3 are determined by 
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solving the following three transcendental equations originating from the boundary condition JSJ: 

[2L -M-N-(i + j)ir] tan L = 3j, < L < 
[2M - N -L + iir] tan A/ = 37, < M < - 

[2N-L-M + jn] tan A = 37, < N < | , (A5) 

where 7 = Rp/D, i = 0, 1, . . ., j = £, £ + 1, . . . and the auxiliary variables L, M and N are related to /1 , z/ , #1,2,3 in the 
following way: 

61 = L — M — N , S 2 = -L + M - N , 6 3 = -L-M + N 

2M-N-L 2N-L-M , , A , 

/i = h i , f = 1- j , A = —fx — v . (A6) 

7T 7T 



Finally, the eigenvalues are given by: 

T.t 1 = 

27 Vi? 



^ 1 = ^^) 2 [A 2 +/x 2 + , 2 ]. (A7) 



N s ' a2 = / / _ dxdyT s , a {x,y)T s , a {x,y) (A8) 

JO Jy/V3 

in an obvious notation, the subscript i,j has been suppressed. After a rather lengthy manipulation, the result can be 
written: 

N-j 2 = N s2 = N a2 = F\jm, v, S u 8 2 ] + F[u, /1, £3] + -/1 - f, <5 2 , <J 3 ] 

+0(A*,*a]+Q[i/,53] + 0[-A*-f.*i] (A9) 



and for i = j: 



A* 2 = ~ 8 - 8/1 V + 7cos[2<5 2 ] + 8cos[2H 

+ cos[2<S 2 - 4/wr] - 8 cos[2<5 2 - 2/«r] - 4/i7r sin[2<5 2 ] + 16^7r sin[25 2 - 2/z7r] } (A10) 
The functions F and Q are given below: 

3\/3-R 2 

F\pi,v, 81,82] = —. — 77 — ; — 7— 5- {^cos[<5i - 6 2 ]0 + /icos[2(/x + i/W]) 
4/xi/ J (// + v)tt z 

+ (p + v)(—v cos[^i - <5 2 + 2/i7r] - ^cos[<5i + 5 2 ] 

+/icos[(5i + <5 2 + 2i/7r] + 2/i^7rsin[<5i + <5 2 ]) — sin[<5i — 8 2 ] sin[2(/j + v)tt]} (All) 
3 /3i? 2 

Qt^. = a o (cos[2(J 2 - M7 r)] - cos[2<5 2 ] + 2 M 7r( M 7r - sin[2(<5 2 - /wr)])) (A12) 

8/i Z 7T Z 

Note that when £ = j we have M = N, S 2 — S 3 , fi — v and 27r/i = <5 2 — <5i + 2i7r. The eigenf unctions and transcendental 
equations simplifies and can after some manipulations be reduced the equations given in section ITTT1 
To summarize the Green's function for the general case can be written : 

- T ^yyr^jx' lV ') 



G(x, y\x', y'\ t) = ^ " ' 
00 00 

+ E E 1^2 [Tt^T^x'^+Tt^x^T^x^y')] (A13) 

i=0 j=i+l ij 

This function can be used to calculate the magnetic signal when there are gradients present in the system. 

In section ITTT1 we used the fact that only the symmetric diagonal modes contribute to the magnetic signal. That 
the antisymmetric modes does not contribute follows by definition, but the non diagonal symmetric modes (£ 7^ j) 
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may give a non-vanishing contribution to the magnetic signal. To show this we integrate the symmetric modes over 
the equilateral triangular domain : 

J-3R r-y^+2V3R 27\/3i? 2 

A= / / dxdyT s (x,y) = -— ——————{(ii-v)cos[S 1 } + ( f i + 2v)cos[S 2 } 

Jo JyV3 27r (P - v) (2/x + (/i + 2i/) 

27T 27T 27T 

-(2/U + i/) cos[(5 3 ] + (\i + 2v) cos[<5 3 + —(//- v)] - (2/i + i/) cos[<5 2 + — (-// + z/)] + (/z - z/) cos[<5 2 — + f)] 

o o o 

27T 27T 27T 

+0 + 2v) cos[5i + —(2^ + i/)] + (/i - z/)cos[<5 3 - — (// + 2^tt)] - (2/i + cos[<5i + — + 2;/)]} (A14) 
o o o 



From equation i|A6(l . we find: 

Si - 2S 2 + S 3 . 
" = 2~, +Z 

<5i — 2<5 3 + 5-) , . ■. 

y= - 27T 3 (A15) 

We need to consider two cases separately, j — i, j = i mod 3(i =/= j). For j = i we have in addition S 3 = S 2 ((J, = v). 
We then find: 



Q\/3/? 2 

A** 1 = -f— -{cosfel - cosfe + 2i7r - 2 M 7r] + 2/z7rsin[<5 2 ]} (A16) 

2[l Z TX jL 



For j = t mod 3(j 7^ i), we write j = i + 3fc, fc = 1, 2, . . 

1SV3R 2 



A (S 2 -6 a - 2kn)(-5 1 + 5 2 + 2(i + fc)7r)(-<5i + <5 3 + 2(* + 2fc)7r) 
x{(S 2 + S3 + 2fc7r) cos[<5i] + (Si -5 3 - 2(i + 2fc)?r) cos[<5 2 ] + (Si + S 2 + 2(i + k)n) cos[c5 3 ]} (A17) 

For the other cases, we find A = 0. Using the constraints imposed by the boundary conditions (|A5J| it turns out that 
A k = 0, this has been verified numerically by solving equation l|A5|l for different values of 7 = pr/D. We are then left 
with equation (|A16|) as the final result. 

APPENDIX B: RANDOM WALK 

1. The relative importance of relaxation types 

Let 9 denote the ratio of bulk relaxation and boundary relaxation rates. In the fast diffusion li mit, with a compara- 
tively flat distribution of walkers, 9 can be estimated by combining the ratio of bulk lattice points (N(N + 1)/2) to 
boundary lattice points (3(-/V— 1)) with the appropriate ratio of relaxation probabilities (a boundary walker will suffer 
relaxation with probability £/3 on this lattice except in the very corner positions, where the probability is 2£/3): 

N 1 + 1/N k_ 1 a N 2 (N + 1) 
T (1 - l/iV)2 C ~ 4y/3pT 2b (i\T-l) 3 

> J^^_1^ZL_ ( Bi) 



4V3pT 2fc 6 7 T, 



2b 



As follows from inspection of equation Ij23(l and (|27(l , and implicitly stated by Brownstein and Tan-HQ 

the solutions 

for the average magnetization M will be three-parameter functions. In addition to 7 and 9, one of the intrinsic time 
scales, td for diffusion and t p for relaxation at the boundary, respectively, may be chosen as a parameter: 

a 2 

TD = W 

( a \2 (2 X ltT 5 cm 2 /s\ 

= 1.25 s X '—) B2) 

V 0.1 mm./ V D j v ' 
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FIG. 9: (N, £)-dependency for 7 = 1, 10 6 walkers (Mersenne Twister). Note that as N increases the numerical solution moves 
towards the analytical solution 




One such set of parameters may thus be 

U '4D' ' 

From equation (IB 1|) it follows that, physically, 6 depends on the bulk sink strength density measured in units of the 
inverse time scale for relaxation at the boundary. 



2. Sources and estimates of simulation inaccuracy 

Figure shows the results of simulations of the average magnetization M for 7 = 1 , in an interval of scaled time 
t/To. The N values 26, 51, 101 and 201 were used; for each, 10 simulations with different random number generator 
seeds have been plotted, with N w = 10 6 random walkers released on the lattice in each simulation. The lower most 
line shows the analytical solution in the same interval, plotted as the sum of 50 modes. The 'Mersenne Twister' 
generator|22j was used in this plot as random number generator, but also runs with single-precision versions of ran2 
and ran3|23j were made. The simulations organize themselves in bands corresponding to the four N values used, with 
N increasing from above (the bands for N = 101 and N = 201 partially merging). 

Various sources of error in the simulations can be discerned: 

• Random: 

— Small fluctuations in the curves increase with increasing N, for the same number of random walkers. 

— Generator seed / number of walkers: With each curve showing an average over N w — 10 6 walkers, for each 
N value there is still a seed-dependent spread within a band of the order of 0.7 %. 

— Truncation errors do not give important random effects in these simulations; curves for a single-precision 
ran3 (not shown) have about the same variance as those shown here. 

• Systematic: 

— Finite e effects: Since e oc 1/(N— 1), choosing N too small makes the order e terms in equation (|25|) increase 
in importance, thus violating equation (|27|l in addition to introducing inaccuracies in the representation of 
the derivatives in equation l|23|) and (|27|l . The simulation results averaged over seed values show a deviation 
from the analytical solution as a monotonous function of N, of order 1 % for N — 101 and rapidly 
increasing as N decreases. (Technically, these iV-dependent systematic deviations may be considered a 
'scaling violation'.) 
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— Truncation errors: For N = 26 the results for ran3 (not shown) and for Mersenne Twister are very close, 
but for increasing N the average difference between the bands becomes larger for ran3 than for Mersenne 
Twister, so that for large N the bands actually make an 'undershoot' (not shown) of the analytical result. 
Such effects may arise from truncation when a random generator is used to make a choice with a given 
probability. 

The level of accuracy may vary with 7 and with the range of scaled time used. For the values treated above, we 
conclude that N ~ 200 and N w ~ 10 6 should be used to obtain an accuracy in one given run of order 1 % or better. 
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